################################################################################
# Created By:Pietryka
# Creation Date:  2016-08-22
# Purpose: Create Table 1, Regressions of Turnout in CPS data
# Questions: mpietryka@fsu.edu
################################################################################


# PREAMBLE  ====================================



# LOAD PACKAGES --------------
library(tidyverse)



# LOAD DATA --------------


# load objects: 'cps_sub_10', 'cps_sub_14', 'main_vars', and 'covariates'
# created in 'CPS-1A-Clean-Turnout.R'
load("Data/CPS-1A-Clean-Turnout.Rdata")




# CLEAN DATA   ====================================

cps_combined <- cps_sub_10  %>%
  tbl_df()  %>%
  mutate(year = 2010)  %>%
  # RECODE RESIDENTIAL DURATION VARIABLE TO MATCH 2014
  mutate(PES8 = recode(PES8,
                       `1` = 1,
                       `2` = 1,
                       `3` = 1,
                       `4` = 2,
                       `5` = 3,
                       `6` = 4))  %>%
  bind_rows(cps_sub_14)  %>%
  # NEW VARIABLES
  mutate(year = ifelse(is.na(year), 2014, year),
         # EMPLOYMENT
         work_employed =  if_else(PEMLR == 1, 1, 0, missing = NA_real_),
         work_retired  =  if_else(PEMLR == 5, 1, 0, missing = NA_real_),
         # RELATIONSHIP STATUS
         marry_married  =  if_else(PEMARITL == 1, 1, 0, missing = NA_real_),
         marry_sepordiv  =  if_else(PEMARITL %in% c(2, 4, 5) ,
                                    1,
                                    0,
                                    missing = NA_real_),
         marry_widow  =  if_else(PEMARITL == 3, 1, 0, missing = NA_real_),
         # EDUCATION
         educ_hs = if_else(educ == 39, 1, 0, missing = NA_real_),
         educ_somecol = if_else(educ > 39 & educ < 43,
                                1,
                                0,
                                missing = NA_real_),
         educ_colgrad = if_else(educ == 43, 1, 0, missing = NA_real_),
         educ_advanced = if_else(educ > 43 & educ <= 46,
                                 1,
                                 0,
                                 missing = NA_real_),
         # RACE
         race_black = if_else(PTDTRACE == 2, 1, 0, missing = NA_real_),
         race_asian = if_else(PTDTRACE == 4, 1, 0, missing = NA_real_),
         race_other = if_else(PTDTRACE %in% c(3, 5, 6:21), 1, 0, missing = NA_real_),
         # RESIDENTIAL MOBILITY
         res_1t2 = if_else(PES8 == 2, 1, 0, missing = NA_real_),
         res_3t4 = if_else(PES8 == 3, 1, 0, missing = NA_real_),
         res_5plus = if_else(PES8 == 4, 1, 0, missing = NA_real_)
  )  %>%
  # FAMILY INCOME
  group_by(year)  %>%
  mutate(famincquart = ntile(famincome, 4))  %>%
  ungroup()


# PREPARE MODELS   ====================================


formula1 <- "turnout ~ prekkids"
formula2 <- paste(formula1,
                  "age",
                  "female",
                  sep = " + ")
formula3 <- paste(formula2,
                  "work_employed",
                  "work_retired",
                  "student",
                  "educ_hs",
                  "educ_somecol",
                  "educ_colgrad",
                  "educ_advanced",
                  "marry_married",
                  "marry_sepordiv",
                  "marry_widow",
                  "factor(famincquart)",
                  "res_1t2",
                  "res_3t4",
                  "res_5plus",
                  "race_black",
                  "race_asian",
                  "race_other",
                  sep = " + ")

year_form <- crossing(year = c(2010, 2014),
                      form = c(formula1, formula2, formula3)
)

fit_logit <- function(the_year, the_formula) {
  data <-   cps_combined  %>% filter(year == the_year)
  glm( formula(the_formula), data =  data, family = "binomial")
}


# FIT MODELS   ====================================
library(texreg)

cps_mods <- apply(year_form, 1,
                  function(x) fit_logit(the_year = x["year"],
                                        the_formula = x["form"])
)






var_map <- list(
  prekkids = "Parent of Young Child (0 = No; 1 = Yes)",
  age = "Age (in Years)",
  female = "Female (0 = No; 1 = Yes)",
  student = "Student (0 = No; 1 = Yes)",
  work_employed = "Employed (0 = No; 1 = Yes)",
  work_retired = "Retired (0 = No; 1 = Yes)",
  educ_hs = "HS Graduate (0 = No; 1 = Yes)",
  educ_somecol = "Some College (0 = No; 1 = Yes)",
  educ_colgrad = "College Graduate (0 = No; 1 = Yes)",
  educ_advanced = "Advanced Degree (0 = No; 1 = Yes)",
  marry_married = "Married (0 = No; 1 = Yes)",
  marry_sepordiv = "Separated/Divorced (0 = No; 1 = Yes)",
  marry_widow = "Widowed (0 = No; 1 = Yes)",
  `factor(famincquart)2` = "In Second Quartile (0 = No; 1 = Yes)",
  `factor(famincquart)3` = "In Third Quartile (0 = No; 1 = Yes)",
  `factor(famincquart)4` = "In Fourth Quartile (0 = No; 1 = Yes)",
  res_1t2 = "1-2 Years (0 = No; 1 = Yes)",
  res_3t4 = "3-4 Years (0 = No; 1 = Yes)",
  res_5plus = "5+ Years (0 = No; 1 = Yes)",
  race_black = "Black (0 = No; 1 = Yes)",
  race_asian = "Asian (0 = No; 1 = Yes)",
  race_other = "Other (0 = No; 1 = Yes)",
 `(Intercept)` = "Intercept"
)

screenreg(cps_mods[c(3, 6)],
          custom.coef.map = var_map)

group_list <- list(
  " " = 1:4,
  "Employment Status (Reference = Unemployed/Not in Labor Force)" = 5:6,
  "Education (Reference = Less than HS Graduate)" = 7:10,
  "Relationship (Reference = Never Married)" = 11:13,
  "Income (Reference = In First Quartile)" = 14:16,
  "Time in Residence (Reference = Less than One Year)" = 17:19,
  "Race (Reference = White)" = 20:22,
  " " = 23
)


screenreg(cps_mods[c(3, 6)],
          custom.model.names = c("2010", "2014"),
          custom.coef.map = var_map,
          groups = group_list,
          ci.force = TRUE)

htmlreg(cps_mods[c(3, 6)] ,
        caption = "Table 1. Regressions of Turnout on Parenthood Status and Controls using CPS Data.",
        custom.model.names = c("2010", "2014"),
        caption.above = TRUE,
        stars = 0.05,
        groups = group_list,
        custom.coef.map = var_map,
        single.row = TRUE,
        custom.note = "* p < 0.05. Source: CPS 2010 and 2014 Voting and Registration Supplements",
        file = "Tables/table-CPS-turnout.html"
)

